---
title: "Replication: Do Voters Reward Value Over Prevention? Evidence from Disaster Policy Preferences"
author:
  - name: Felix Hartmann
    orcid: 0000-0002-6134-7126
    email: feha.egb@cbs.dk
    affiliations:
      - name: Copenhagen Business School
date: "`r format(Sys.time(), '%d %B, %Y')`"
format:
  html:
    toc: true
    toc-location: left
    code-fold: show
    code-tools: true
    embed-resources: true
---

# Summary information (README)

```{r class.source = "fold-show", eval = FALSE, echo = TRUE}
----------------------------
initial deposit date: 08/01/2025
-----------------------------
Replication code and data needed to reproduce results presented in-text and in the Supplemental Information provided for the study. 

-----------------------------
# code overview
- "0_replication.qmd" R Quarto script produces all the numeric results, tables, and figures in the main text and in the supplementary information

-----------------------------
# files:
- "1_input" contains all data needed for replication
- "2_output" contains numeric results, tables, and figures 

-----------------------------
# 1_input:
- "1_input/afrobarometer_release-dataset_merge-34ctry_r8_en_2023-03-01.sav"
  + Afrobarometer Merged Round 8 data (34 countries) (2022)
  + source: https://www.afrobarometer.org/data/merged-data/

- "1_input/emdat_public_2020_12_11_query_uid-GZ254X.xlsx"
  + EM-DAT: The CRED/OFDA International Disaster Database
  + source: https://public.emdat.be/ 
  
- "1_input/ND_Gain_id_infr_06_raw0.csv"
  + Notre Dame Global Adaptation Initiative Country Index (ND-GAIN)
  + variable: ID_INFR_06 (disaster preparedness)
  + source: https://gain.nd.edu/our-work/country-index/download-data/
  
- "1_input/custom.geo.json"
  + Africa shapefile

- "1_input/conjoint.rds" 
  + contains all data from the conjoint experiment 
	+ long format 
- "1_input/covariates.rds"
  + contains all covariates  
  + merge by case id  to conjoint data
- "1_input/euspeech.korpus.RData" 
  + EUSpeech dataset (Schumacher et al., 2016)
- "1_input/policy_agendas_english.RData" 
  + source: https://github.com/cbpuschmann/quanteda_mzes/tree/master 


-----------------------------
# running time 
- For "0_replication.qmd", about 5 min in total

```
# Housekeeping

## load packages
```{r setup, echo=TRUE}
rm(list=ls())  

options(modelsummary_format_numeric_latex = "mathmode")

knitr::opts_chunk$set(echo = TRUE)

if (!require(pacman)) install.packages("pacman")
pacman::p_load(
  car,   # linear hypothesis check 
  dplyr,
  tidyverse,
  estimatr,
  here,
  reshape,
  quanteda, # speech data
  readtext, # speech data
  lubridate,# speech data
  cregg, # subgroup analysis 
  compare,
  modelsummary,
  texreg,
  gridExtra,
  grid,
  ggplot2,
  lattice,
  vtable,   # Summary Stats
  broom,   # Summary Stats
  broom.helpers,
  lubridate,
  fBasics,
  readr,
  ggpubr,
  foreign,
  readxl, # read excel files
  viridisLite, # colour scheme
  viridis,# colour scheme
  sf, #  spatial data
  ggthemes 
)

```

## load function
```{r}
######################################
# functions
######################################
tidy_coefs_with_ref <- function(mod_obj, sep = "_"){

  tidy_coefs <- tidy(mod_obj) %>% 
    separate(term, c("variable", "level"), sep, remove = FALSE) %>% 
    mutate(variable = paste0(variable, sep))

  xlevels <- mod_obj$xlevels  

  missing_levels <- xlevels %>% 
    enframe() %>% 
    unnest(cols = c(value)) %>% 
    set_names(c("variable", "level"))

  missing_levels %>% 
    anti_join(tidy_coefs) %>% 
    bind_rows(tidy_coefs) 
}
```



## Session Info
```{r}
sessionInfo()

```
## load data
```{r}
here::i_am("0_replication.qmd")
df_long <- readRDS("1_input/conjoint.rds")
df_long<-as_tibble(df_long)
```


# Main Analysis
## Figure 1: Main Results
```{r,warning=FALSE}
df_long_main<-df_long %>% dplyr::rename(
                          EffortPrevention_ = Effort_prevention,
                          EffortRelief_ = Effort_relief,
                          EffectivePrevention_ = Effective_prevention,
                          EffectiveRelief_ = Effective_relief,
                          Ask_ = Ask,
                          EQ_ = EQ,
                          Honesty_ = Honesty
                          )


  res<-lm_robust(Choice ~EffortPrevention_+EffortRelief_+EffectiveRelief_+EffectivePrevention_+Ask_+EQ_+Honesty_,
                  data=df_long_main,clusters = CaseID, se_type = "stata") %>%
  tidy_coefs_with_ref() %>%
  # remove intercept
  dplyr::filter(!(variable == '(Intercept)_')) %>%
  # add space for level
  dplyr::mutate(level = paste(level, "   ")) %>%
  # fill value
  dplyr::mutate(estimate = tidyr::replace_na(estimate, 0)) %>% 
    # rename
  dplyr::mutate(	
  variable=dplyr::recode(variable, 
                     "Ask_"="E. Ask_", 
                     "Honesty_" = "G. Honesty_",
                     "EQ_" = "F. EQ_",
                     "EffectivePrevention_" = "C. EffectivePrevention_",
                     "EffectiveRelief_" = "D. EffectiveRelief_",
                     "EffortPrevention_" = "A. EffortPrevention_",
                     "EffortRelief_" = "B. EffortRelief_"
  ))    %>% 
  # add empty raw
  as_tibble() %>%
  group_by(variable) %>% 
  group_modify(~ add_row(.x,.before=0)) %>%
  #
  dplyr::mutate(level = fct_reorder(level, estimate,.na_rm = TRUE))  %>% 
  dplyr::mutate(level = dplyr::coalesce(level,variable))   %>% 
    # rename
  dplyr::mutate(	
  level=dplyr::recode(level, 
                     "A. EffortPrevention_" = "Effort Preparedness",
                     "B. EffortRelief_" = "Effort Relief",
                     "C. EffectivePrevention_" = "Effective Preparedness",
                     "D. EffectiveRelief_" = "Effective Relief",
                     "E. Ask_"="Ask Help", 
                     "F. EQ_" = "Visits",
                     "G. Honesty_" = "Corruption"
  )) %>%
  add_column(indicator = 0) %>%
  # replace
  dplyr::mutate(
    indicator=replace(indicator, variable=="A. EffortPrevention_", "Effort"),
    indicator=replace(indicator, variable=="B. EffortRelief_", "Effort"),
    indicator=replace(indicator, variable=="C. EffectivePrevention_", "Effective"),
    indicator=replace(indicator, variable=="D. EffectiveRelief_", "Effective"),
    indicator=replace(indicator, variable=="E. Ask_", "Other"),
    indicator=replace(indicator, variable=="F. EQ_", "Other"),
        indicator=replace(indicator, variable=="G. Honesty_", "Other")
      )  %>%
  dplyr::mutate(	
  level=dplyr::recode(level, 
                     "Prevention Effective    "="Yes    ", 
                     "Prevention Not Effective    " = "No    ",
                     "Relief Effective    " = " Yes    ",
                     "No Relief    " = " No    ",
                     "Prevention Effort    " = "High    ",
                     "No Prevention Effort    " = "Low    ",
                     "Relief Effort    " = " High   ",
                     "No Relief Effort    " = " Low    ",
                     "Ask Relief    " = "  Yes    ",
                     "Not Ask Relief    " = "  No    ",
                     "No Visits    " = "   No    ",
                     "Visits    " = "   Yes    "

  ))



res$level <- factor(res$level, levels = res$level)
res<- res %>%  mutate(level = fct_rev(as_factor(level))) 


bold.ticks <- c("Ask Help","Corruption","Visits","Effective Preparedness","Effective Relief","Effort Preparedness","Effort Relief")
bold.labels <- ifelse(levels(res$level) %in% bold.ticks, yes = "bold", no = "plain")

fig1<-res %>%  
  ggplot( aes(x = estimate, y =  level,color = indicator)) +
    geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = level, 
                     xmin = conf.low, 
                     xmax = conf.high),
                 size=0.5, alpha = 1, height = 0.2) +
  scale_colour_colorblind()+
  #scale_fill_viridis_b(name = 'indicator') + 
  theme_bw() +
  ggtitle("AMCEs")+
  xlab("Effect on Probability of Support for Candidate") +
  ylab("Disaster Policy Choice")+
  theme(axis.text=element_text(size = 12)) +
  theme(strip.text.x = element_text(size =12))+
  theme(axis.text.y = element_text(face=bold.labels))+
  theme(axis.title = element_text(size=12,face="bold",vjust=-2.5))+
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")+
  theme(legend.position = "right")

fig1
ggsave("2_output/fig_1.pdf", width = 8, height = 4, units = "in")
```
## Figure 2: Marginal Effect of Prevention Effort Conditional on Effectiveness.
```{r,warning=FALSE}

# Q: Does support for prevention efforts gets lower if voters see it as not being effective?
# create indicator for occurrence of combination across rows

df_long<-  df_long %>%
  dplyr::mutate(
    Effort_prevention_num = case_when(
      Effort_prevention == "No Prevention Effort"  ~ 0,
      Effort_prevention == "Prevention Effort"  ~ 1
    ),
    Effort_relief_num = case_when(
      Effort_relief == "No Relief Effort"  ~ 0,
      Effort_relief == "Relief Effort"  ~ 1
      ),
    Effective_prevention_num = case_when(
      Effective_prevention == "Prevention Not Effective"  ~ 0,
      Effective_prevention == "Prevention Effective"  ~ 1
  ),
    Effective_relief_num = case_when(
      Effective_relief == "No Relief"  ~ 0,
      Effective_relief == "Relief Effective"  ~ 1
  ),
    Ask_num = case_when(
      Ask == "Not Ask Relief"  ~ 0,
      Ask == "Ask Relief"  ~ 1
  ),
    EQ_num = case_when(
      EQ == "No Visits"  ~ 0,
      EQ == "Visits"  ~ 1
  ),
  Corruption_num = case_when(
    Honesty == "No Corruption"  ~ 0,
    Honesty == "Corruption"  ~ 1,
    Honesty == "Vote Buying"  ~ 2
  )
  )

  # count prevention effort + not effective by respondent
  df_long<-df_long %>% 
  dplyr::mutate(prevention_ineffective = dplyr::case_when(Effort_prevention_num == 1 & Effective_prevention_num == 0 ~ 1)) %>% 
  dplyr::mutate(prevention_ineffective = replace_na(prevention_ineffective, 0)) %>% 
  group_by(CaseID) %>% 
  mutate(Frequency_prevention_ineffective = sum(prevention_ineffective))

  df_long$Frequency_prevention_ineffective <- as.factor(df_long$Frequency_prevention_ineffective) 
  df_long$Effort_prevention_num <- as.factor(df_long$Effort_prevention_num) 

  prevention_effort_ineffective<-df_long %>%
  cj(Choice ~ Effort_prevention_num ,id = ~ CaseID, estimate = "amce",by = ~Frequency_prevention_ineffective) %>%  
  dplyr::filter(level == 1)

  prev_ineff<-prevention_effort_ineffective %>%
  ggplot( aes(x = estimate, y = BY)) +
    geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = BY, 
                     xmin = estimate - 1.96*std.error, 
                     xmax = estimate + 1.96*std.error),
                 linewidth=0.5, alpha = 0.5, height = 0.2) +
  scale_colour_viridis_d(name = 'Group') + 
  coord_flip()+
  theme_bw() +
  ggtitle("AMCEs: Prevention Effort")+
  scale_x_continuous(limits = c(-0.5, 1))+
  xlab("Estimate") + ylab("Number of Profiles: High Effort + No Outcome")+
  theme(axis.text=element_text(size = 12)) +
  #theme(axis.text.x = element_text(size = 12)) +
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")


  # Prevention effort + effective
  df_long<-df_long %>% 
  dplyr::mutate(prevention_effective = dplyr::case_when(Effort_prevention_num == 1 & Effective_prevention_num == 1 ~ 1)) %>% 
  dplyr::mutate(prevention_effective = replace_na(prevention_effective, 0)) %>% 
  group_by(CaseID) %>% 
  mutate(Frequency_prevention_effective = sum(prevention_effective))

  df_long$Frequency_prevention_effective <- as.factor(df_long$Frequency_prevention_effective) 
  df_long$Effort_prevention_num <- as.factor(df_long$Effort_prevention_num) 

  
  prevention_eff<-df_long %>% 
  # could exclude contest 8 (too few cases) and cluster by CaseID, results do not change
    cj( Choice ~ Effort_prevention_num ,estimate = "amce",by = ~Frequency_prevention_effective) %>% 
    dplyr::filter(level == 1)

  prev_eff<-prevention_eff %>%
  #mutate(level = fct_reorder(estimate, desc(BY))) %>%
  ggplot( aes(x = estimate, y = BY)) +
    geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = BY, 
                     xmin = estimate - 1.96*std.error, 
                     xmax = estimate + 1.96*std.error),
                 size=0.5, alpha = 0.5, height = 0.2) +
    scale_colour_viridis_d(name = 'Group') + 
  coord_flip()+
  theme_bw() +
  ggtitle("AMCEs: Prevention Effort")+
  scale_x_continuous(limits = c(-0.5, 1))+
  theme(axis.text=element_text(size = 12)) +
  xlab("Estimate") + ylab("Number of Profiles: High Effort + Effective Outcome")+
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")

fig_updating  <- ggarrange(prev_ineff,prev_eff,
          ncol = 2, nrow = 1, common.legend = TRUE, legend="bottom")
fig_updating

ggsave("2_output/fig_2_updating.pdf", width = 9, height = 3, units = "in")
```

# Online Appendix

## Figure A1: Natural Disasters across Africa 1970-2020
```{r,warning=FALSE}
################################################
# Total Number of Disasters by country
################################################

emdat <- read_excel("1_input/emdat_public_2020_12_11_query_uid-GZ254X.xlsx")

# change name 
names(emdat)[names(emdat) == 'Disaster Type'] <- 'Disaster'

# keep most common disaster (scale gets to messy otherwise)
emdat <- emdat[ which(emdat$Disaster=='Flood'  | emdat$Disaster=='Drought' | emdat$Disaster=='Storm' | emdat$Disaster=='Wildfire' | emdat$Disaster=='Insect infestation'| emdat$Disaster=='Extreme temperature'), ]

# subset, 1970 +
emdat<-emdat %>% dplyr::filter(as.numeric(Year) > 1969)
  
# create event variable for disaster
emdat$total <- rep(1,nrow(emdat)) 

# aggregate by group year 
DisasterAggTypeYear<-emdat %>% 
  group_by(Year, Disaster) %>% 
  summarise_at(vars(total), list(~sum(.)))

# Plot
fig_disaster_time<-ggplot() + geom_bar(aes(y = total, x =as.numeric(Year), fill = Disaster),
                      data = DisasterAggTypeYear,
                      stat="identity") + 
  scale_x_continuous(breaks = seq(1970, 2020, by = 10)) +
  theme_bw(base_size = 12) +
  scale_fill_viridis_d() +
  #scale_fill_grey()+
  theme(legend.position="bottom",
        legend.justification="bottom",
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(-10,-10,-10,-10))+
  labs(x = "Year", y ="Number of Disasters")+
  ggtitle("Africa: Disasters")
  
##################################
# Malawi data
##################################

Malawi <- emdat[ which(emdat$Country=='Malawi'), ]
Malawi_flood <- Malawi[ which(Malawi$Disaster=='Flood'), ]

# create event variable for disaster
emdat$total <- rep(1,nrow(emdat)) 

Malawi_agg_flood<-Malawi_flood %>% 
  group_by(Year) %>% 
  summarise_at(vars(total), list(~sum(.)))

Malawi_agg_year<-Malawi %>% 
  group_by(Year) %>% 
  summarise_at(vars(total), list(~sum(.)))

# fill time series gaps
Year<-seq(1960, 2020, 1,)
Full_years<- as.data.frame(Year)
Malawi_agg_year$Year <- as.numeric(Malawi_agg_year$Year)
Malawi_agg_flood$Year <- as.numeric(Malawi_agg_flood$Year)

# merge
Malawi_agg_year_full<- Full_years %>% left_join(Malawi_agg_year,by=c("Year"))
Malawi_agg_flood_full<- Full_years %>% left_join(Malawi_agg_flood,by=c("Year"))

# fill in NA as 0
Malawi_agg_year_full <- Malawi_agg_year_full %>%
  mutate(total_full = if_else(is.na(total), 0, total))

Malawi_agg_flood_full <- Malawi_agg_flood_full %>%
  mutate(flood = if_else(is.na(total), 0, total))


Malawi_final<- Malawi_agg_year_full %>% left_join(Malawi_agg_flood_full,by=c("Year"))


# create three year moving averages, otherwise not so informative 
Malawi_final <- Malawi_final %>%
  dplyr::select(Year, srate = total_full, srate2 = flood) %>%
  dplyr::mutate(
    All = zoo::rollmean(srate, k = 5, fill = NA, align = "center"),
    Flood = zoo::rollmean(srate2, k = 5, fill = NA, align = "center")
  )

# melt 
Malawi_final<-melt(data = Malawi_final, id.vars = "Year", measure.vars = c("All", "Flood"))

# change names
names(Malawi_final)[names(Malawi_final) == "variable"] <- "Disaster"


#Plot by year
fig_disaster_malawi_time<-ggplot(data = Malawi_final) +
  geom_line(aes(y = value, x = Year , colour=Disaster), 
            alpha = 0.6,
            size = 0.6) +
  scale_fill_viridis() +
  #scale_color_manual(values = c("gray40","#09557f")) +
  theme_bw(base_size = 12) +
  theme(legend.position="bottom",
        legend.justification="bottom",
        legend.margin=margin(0,0,0,0),
        legend.box.margin=margin(-10,-10,-10,-10))+
  scale_x_continuous(breaks = seq(1970, 2020, by = 10)) +
  scale_y_continuous(breaks = seq(0, 5)) +
  labs(x = "Year", y ="Number of Disasters")+
  ggtitle("Malawi: Disasters")
  #xlim(c(1965,2020))+ 
  #scale_fill_discrete(name = "New Legend Title")
  
fig_disaster_timeseries<-grid.arrange(fig_disaster_time,fig_disaster_malawi_time , ncol = 2)

fig_disaster_timeseries
ggsave("2_output/fig_A1_disaster_timeseries.pdf",width = 10, height = 3, units = "in",fig_disaster_timeseries )


```

## Figure A2: Number of Natural Disasters by Country (1970 and 2020) and Preparedness by Country (2007-2011).

```{r}

# Read GeoJSON file
africa_shp <- st_read("1_input/custom.geo.json")
# Extract iso_a3 values from africa_shp
iso_a3_values <- africa_shp %>% 
  st_drop_geometry() %>%  # Remove geometry to just keep the data
  pull(iso_a3)  # Extract iso_a3 column


# path in downloaded data: ND_Gain/resources/indicators/id_infr_06/raw0.csv
africa_data <- read.csv(file = '1_input/ND_Gain_id_infr_06_raw0.csv') %>%
  mutate(mean_col = rowMeans(select(., X2007, X2009, X2011), na.rm = TRUE)) %>%
  select(ISO3, Name, X2007, X2009, X2011, mean_col)%>% 
  filter(ISO3 %in% iso_a3_values)

# calculate points at which to plot labels
centroids <- africa_shp %>% 
  st_centroid() %>% 
  bind_cols(as_tibble(st_coordinates(.)))    # unpack points to lat/lon columns

centroids <- centroids %>%
  mutate(country.flag = ifelse(iso_a3 == "MWI", "MW", NA)) # "MW" for Malawi

spatial_preparedness_africa<-africa_data %>% 
  left_join(africa_shp, ., by = c('iso_a3' = 'ISO3')) %>% 
  ggplot() + 
  geom_sf(aes(fill = mean_col)) + 
  scale_fill_viridis( name="Disaster Preparedness" ) +
  geom_text(aes(X, Y, label = country.flag), data = centroids, size = 4, color = "red", na.rm = TRUE) +
  theme_bw(base_size = 12) +
  xlab("Latitude") + 
  ylab("Longitude")+
  #ggtitle("Africa")+
  theme(legend.position="bottom")



################################################
# Total Number of Disasters by country
################################################


agg_country<-emdat %>% 
  group_by(ISO) %>% 
  summarise_at(vars(total), list(~sum(.)))

spatial_conflict_africa<-agg_country %>% 
  left_join(africa_shp, ., by = c('iso_a3' = 'ISO')) %>% 
  ggplot() + 
  geom_sf(aes(fill = total)) + 
  scale_fill_viridis(name="Total Disasters" ) +
  geom_text(aes(X, Y, label = country.flag), data = centroids, size = 4, color = 'red')+
  theme_bw(base_size = 12) +
  xlab("Latitude") + 
  ylab("Longitude")+
  #ggtitle("Africa: Total Number of Natural Disasters")+
  theme(legend.position="bottom")



fig_A2_spatial<-grid.arrange(spatial_conflict_africa,                             
                    spatial_preparedness_africa,
                    nrow = 1,ncol=2)        
ggsave("2_output/fig_A2_spatial.pdf",width = 10, height = 5, units = "in",fig_A2_spatial )


```



## Figure A3: Malawi Disaster Spending and Perceived Responsibilities
```{r,warning=FALSE}

################################################
################################################
# Relief Aid 
################################################
################################################

# Humanitarian Responses: https://www.usaid.gov/sites/default/files/documents/1860/Malawi_National_Resilience_Strategy.pdf
# Page: 10
Relief_million <- c(2.5, 
                    5, 
                    7,
                    25,
                    20,
                    27.5,
                    55,
                    40,
                    57.5,
                    57.5,
                    57.5,
                    62,
                    45,
                    30,
                    80,
                    100,
                    350)

Year<-seq(from = 2000, to = 2016, by=1)


response <- as.data.frame(cbind(Year,Relief_million))


##########################################
# Adaption (Prevention id)
##########################################

Malawi_Geocoded_Climate_Dataset <- read_excel("1_input/Malawi Geocoded Climate Dataset.xlsx", 
                                              sheet = "All Projects")

climate_projects <- subset(Malawi_Geocoded_Climate_Dataset, `Climate Spectrum Assignment` =="Climate Oriented")

# keep if amount specified
climate_projects<-climate_projects[!is.na(climate_projects$`Cumulative Disbursement`), ]

# if date of signment is "NA" or "0", use the completion date if possible
climate_projects$`Date of Agreement Signed` <- ifelse(is.na(climate_projects$`Date of Agreement Signed`), climate_projects$`Date of Planned Completion`, climate_projects$`Date of Agreement Signed`)
climate_projects$`Date of Agreement Signed` <- ifelse(climate_projects$`Date of Agreement Signed`==0, climate_projects$`Date of Planned Completion`, climate_projects$`Date of Agreement Signed`)

# keep if year of agreement is identified 
climate_projects<-climate_projects[!is.na(climate_projects$`Date of Agreement Signed`), ]

# drop if year is marked as "0" (coding error?)
climate_projects<-subset(climate_projects,`Date of Agreement Signed`!=0)

# amount in millins US $
climate_projects[c("Cumulative Disbursement")] <- climate_projects[c("Cumulative Disbursement")]/1e6


# standardize date
climate_projects<-climate_projects %>% 
  mutate(Year = as.character(`Date of Agreement Signed`)) %>% 
  mutate(Year = replace(Year, Year == '01/03/2006', '2006')) %>% 
  mutate(Year = replace(Year, Year == '01/05/2009', '2009')) %>% 
  mutate(Year = replace(Year, Year == '01/07/1999', '1999')) %>% 
  mutate(Year = replace(Year, Year == '01/07/2001', '2001')) %>% 
  mutate(Year = replace(Year, Year == '01/07/2002', '2002')) %>% 
  mutate(Year = replace(Year, Year == '01/07/2003', '2003')) %>% 
  mutate(Year = replace(Year, Year == '01/07/2006', '2006')) %>% 
  mutate(Year = replace(Year, Year == '01/07/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '01/07/2010', '2010')) %>% 
  mutate(Year = replace(Year, Year == '01/08/2005', '2005')) %>% 
  mutate(Year = replace(Year, Year == '01/09/2007', '2007')) %>% 
  mutate(Year = replace(Year, Year == '01/09/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '01/10/2004', '2004')) %>% 
  mutate(Year = replace(Year, Year == '01/11/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '01/11/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '01/12/2007', '2007')) %>% 
  mutate(Year = replace(Year, Year == '01/12/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '02/01/2001', '2001')) %>% 
  mutate(Year = replace(Year, Year == '02/01/2003', '2003')) %>% 
  mutate(Year = replace(Year, Year == '11/04/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '12/09/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '12/09/2008', '2008')) %>% 
  mutate(Year = replace(Year, Year == '15/12/2006', '2006')) %>% 
  mutate(Year = replace(Year, Year == '23/01/2007', '2007')) %>% 
  mutate(Year = replace(Year, Year == '19/11/2008', '2008'))%>% 
  mutate(Year = replace(Year, Year == '05/05/2010', '2010')) %>% 
  mutate(Year = replace(Year, Year == '31/12/2009', '2009')) %>% 
  mutate(Year = replace(Year, Year == '31/03/2011', '2011')) %>% 
  mutate(Year = replace(Year, Year == '31/03/2010', '2010')) %>% 
  mutate(Year = replace(Year, Year == '30/08/2013', '2013')) %>% 
  mutate(Year = replace(Year, Year == '30/06/2009', '2009'))



# only keep unique
# some monyy targets to many places
# Remove duplicates based on Sepal.Width columns
climate_projects<-climate_projects[!duplicated(climate_projects$`Cumulative Disbursement`), ]


climate_projects<-climate_projects%>% 
  group_by(`Year`) %>% 
  summarise_at(vars(`Cumulative Disbursement` ), funs(sum))

# Merge Prevention and Relief Aid Spending 

# INNER JOIN: returns rows when there is a match in both tables.
prevention_response<-merge(climate_projects, response, all.y=TRUE)

# rename
names(prevention_response)[names(prevention_response) == "Cumulative Disbursement"] <- "Prevention"
names(prevention_response)[names(prevention_response) == "Relief_million"] <- "Relief"

# melt 
prevention_response_final<-melt(data = prevention_response, id.vars = "Year", measure.vars = c("Prevention", "Relief"))


# Plot
fig_spending<-ggplot() + geom_bar(aes(y = value, x =as.numeric(Year), fill = variable),
                        data = prevention_response_final,
                        stat="identity")+ 
  theme_bw(base_size = 16) +
  scale_fill_grey()+
  theme(legend.position="bottom")+
  labs(x = "Year", y ="Total Amount in Million US$")+
  ggtitle("Prevention & Relief Spending")+
  guides(fill=guide_legend(title=" "))




# Responsabilities
covariates <- readRDS("1_input/covariates.rds")
covariates <- covariates %>%
  dplyr::rename(
    "resp_prepar_ta" = T_q38_1,
    "resp_prepar_village_head" = T_q38_2,
    "resp_prepar_relig_leader" = T_q38_3,
    "resp_prepar_mp" = T_q38_4,
    "resp_prepar_local_council" = T_q38_5,
    "resp_prepar_district_commissioner" = T_q38_7,
    "resp_prepar_ngo_io" = T_q38_8,
    "resp_prepar_dodma" = T_q38_9,
    "resp_relief_ta" = T_q37_1,
    "resp_relief_village_head" = T_q37_2,
    "resp_relief_relig_leader" = T_q37_3,
    "resp_relief_mp" = T_q37_4,
    "resp_relief_local_council" = T_q37_5,
    "resp_relief_district_commissioner" = T_q37_7,
    "resp_relief_ngo_io" = T_q37_8,
    "resp_relief_dodma" = T_q37_9
  )

# List of the renamed variables for preparation and relief
prep_vars <- c("resp_prepar_ta", "resp_prepar_village_head", "resp_prepar_relig_leader", 
               "resp_prepar_mp", "resp_prepar_local_council", "resp_prepar_district_commissioner", 
               "resp_prepar_ngo_io", "resp_prepar_dodma")

relief_vars <- c("resp_relief_ta", "resp_relief_village_head", "resp_relief_relig_leader", 
                 "resp_relief_mp", "resp_relief_local_council", "resp_relief_district_commissioner", 
                 "resp_relief_ngo_io", "resp_relief_dodma")

# Function to count the number of "3" responses for a given set of variables
count_responsibility <- function(df, variables) {
  # Ensure that the variables are numeric and count the "3" responses
  sapply(variables, function(var) sum(as.numeric(df[[var]]) == 3, na.rm = TRUE))
}


# Count "3" responses for preparation and relief
Prevention <- count_responsibility(covariates, prep_vars)
Relief <- count_responsibility(covariates, relief_vars)

# Create a dataframe for plotting
Office <- c("TA", "Village \n Head", "Relig. \n Leader", "MP", "Council", 
            "District", "IO \n NGO", "DoDMA")
# Combine the counts into a single dataframe
response <- data.frame(Office, Prevention, Relief)

# Reshape the data for plotting (long format)
response_long <- response %>%
  pivot_longer(cols = c("Prevention", "Relief"), names_to = "Category", values_to = "Count")

# Plot the data
responsability<-ggplot(response_long, aes(x = reorder(Office, -Count), y = Count, fill = Category)) +
  geom_col(position = "dodge") +
  scale_fill_grey()+
  theme_bw(base_size = 16) +
  theme(legend.position="bottom")+
  labs(x = "Year", y ="Total Number of Responses")+
  ggtitle("Prevention & Relief Responsibilities Malawi")+
  guides(fill=guide_legend(title=" "))



fig_A3<-ggarrange(
  fig_spending, responsability,
  common.legend = TRUE, legend = "bottom"
  )
fig_A3

ggsave("2_output/fig_A3_malawi.pdf", width = 12, height = 4, units = "in")


```

## Figure A5: Help received after the 2015 flood (left) and satisfaction with response (right)

```{r , echo=TRUE,warning=FALSE}
covariates<-readRDS(file = "1_input/covariates.rds")


covariates<-covariates %>% 
    drop_na(satisfied_personal_help) %>%
     dplyr::mutate(
       personal_help = dplyr::case_when(
              personal_help == 1 ~ "No",
              personal_help == 2 ~ "Yes"
              ),
       community_help = dplyr::case_when(
              community_help == 1 ~ "No",
              community_help == 2 ~ "Yes"
              ),
       satisfied_personal_help = dplyr::case_when(
              satisfied_personal_help == 1 ~ "Very Satisfied",
              satisfied_personal_help == 2 ~ " Somewhat Satisfied",
              satisfied_personal_help == 3 ~ "  Somewhat Dissatisfied",
              satisfied_personal_help == 4 ~ "   Very Dissatisfied",
                         TRUE ~ satisfied_personal_help),
      community_help_id1 = dplyr::case_when(
              community_help_id1 == 1 ~ "Traditional authority",
              community_help_id1 == 2 ~ "Village head",
              community_help_id1 == 3 ~ "Religious leader",
              community_help_id1 == 4 ~ "MP",
              community_help_id1 == 5 ~ "Local council member",
              community_help_id1 == 6 ~ "District commissioner",
              community_help_id1 == 7 ~ "NGO, International Organization",
              community_help_id1 == 8 ~ "Dodma",
              TRUE ~ community_help_id1),
      community_help_id2 = dplyr::case_when(
              community_help_id2 == 1 ~ "Traditional authority",
              community_help_id2 == 2 ~ "Village head",
              community_help_id2 == 3 ~ "Religious leader",
              community_help_id2 == 4 ~ "MP",
              community_help_id2 == 5 ~ "Local council member",
              community_help_id2 == 6 ~ "District commissioner",
              community_help_id2 == 7 ~ "NGO, International Organization",
              community_help_id2 == 8 ~ "Dodma",
              TRUE ~ community_help_id2),
      community_help_id3 = dplyr::case_when(
              community_help_id3 == 1 ~ "Traditional authority",
              community_help_id3 == 2 ~ "Village head",
              community_help_id3 == 3 ~ "Religious leader",
              community_help_id3 == 4 ~ "MP",
              community_help_id3 == 5 ~ "Local council member",
              community_help_id3 == 6 ~ "District commissioner",
              community_help_id3 == 7 ~ "NGO, International Organization",
              community_help_id3 == 8 ~ "Dodma",
                TRUE ~ community_help_id3),
      personal_help_id1 = dplyr::case_when(
              personal_help_id1 == 1 ~ "Family",
              personal_help_id1 == 2 ~ "Friends",
              personal_help_id1 == 3 ~ "Traditional Authority",
              personal_help_id1 == 4 ~ "Village Head",
              personal_help_id1 == 5 ~ "Religious Leader",
              personal_help_id1 == 6 ~ "MP",
              personal_help_id1 == 7 ~ "Local Council",
              personal_help_id1 == 8 ~ "District Commissioner",
              personal_help_id1 == 9 ~ "NGO, International Organization",
              personal_help_id1 == 10 ~ "Department of Disaster Management",
              personal_help_id1 == 11 ~ "Other",
                TRUE ~ personal_help_id1),
      personal_help_id2 = dplyr::case_when(
              personal_help_id2 == 1 ~ "Family",
              personal_help_id2 == 2 ~ "Friends",
              personal_help_id2 == 3 ~ "Traditional Authority",
              personal_help_id2 == 4 ~ "Village Head",
              personal_help_id2 == 5 ~ "Religious Leader",
              personal_help_id2 == 6 ~ "MP",
              personal_help_id2 == 7 ~ "Local Council",
              personal_help_id2 == 8 ~ "District Commissioner",
              personal_help_id2 == 9 ~ "NGO, International Organization",
              personal_help_id2 == 10 ~ "Department of Disaster Management",
              personal_help_id2 == 11 ~ "Other",
                TRUE ~ personal_help_id2),
      personal_help_id3 = dplyr::case_when(
              personal_help_id3 == 1 ~ "Family",
              personal_help_id3 == 2 ~ "Friends",
              personal_help_id3 == 3 ~ "Traditional Authority",
              personal_help_id3 == 4 ~ "Village Head",
              personal_help_id3 == 5 ~ "Religious Leader",
              personal_help_id3 == 6 ~ "MP",
              personal_help_id3 == 7 ~ "Local Council",
              personal_help_id3 == 8 ~ "District Commissioner",
              personal_help_id3 == 9 ~ "NGO, International Organization",
              personal_help_id3 == 10 ~ "Department of Disaster Management",
              personal_help_id3 == 11 ~ "Other",
                TRUE ~ personal_help_id3),
      personal_help_id4 = dplyr::case_when(
              personal_help_id4 == 1 ~ "Family",
              personal_help_id4 == 2 ~ "Friends",
              personal_help_id4 == 3 ~ "Traditional Authority",
              personal_help_id4 == 4 ~ "Village Head",
              personal_help_id4 == 5 ~ "Religious Leader",
              personal_help_id4 == 6 ~ "MP",
              personal_help_id4 == 7 ~ "Local Council",
              personal_help_id4 == 8 ~ "District Commissioner",
              personal_help_id4 == 9 ~ "NGO, International Organization",
              personal_help_id4 == 10 ~ "Department of Disaster Management",
              personal_help_id4 == 11 ~ "Other",
                TRUE ~ personal_help_id4)
      ) 
  

personal_help<-covariates %>%
  dplyr::count(personal_help) %>%
  dplyr::rename(sum_group = n) %>%
  ggplot() +
  geom_bar(aes(x=personal_help,y = (sum_group)/sum(sum_group)),stat="identity",width=0.4) +
  scale_y_continuous(labels=scales::percent) +
  labs(y='Density',x='After the 2015 floods, did you personally receive help?')+
  theme_bw()


satisfied_personal_help<-covariates %>%
  filter(satisfied_personal_help!=".") %>%
  dplyr::count(satisfied_personal_help) %>%
  dplyr::rename(sum_group = n) %>%
  ggplot() +
  geom_bar(aes(x=satisfied_personal_help,y = (sum_group)/sum(sum_group)),stat="identity",width=0.4) +
  scale_y_continuous(labels=scales::percent) +
  labs(y='Density',x='How satisfied were you with the help you received?')+
  theme_bw()



help_rec <- ggarrange(personal_help,satisfied_personal_help,
                    ncol = 2, nrow = 1)
help_rec

ggsave("2_output/fig_A5_help_rec.pdf", width = 12, height = 4, units = "in")


```

## Figure A6: Actors help received (left) and satisfaction when received help only from MP (right)

```{r , echo=TRUE,warning=FALSE}

################################################
################################################
# [q31/Q_48_S] From whom did you receive help? (only personal)
################################################
################################################
personal_help_id1<-as.data.frame(covariates$personal_help_id1)
personal_help_id2<-as.data.frame(covariates$personal_help_id2)
personal_help_id3<-as.data.frame(covariates$personal_help_id3)
personal_help_id4<-as.data.frame(covariates$personal_help_id4)

# drop all missing
personal_help_id1<-subset(personal_help_id1, personal_help_id1!="." & personal_help_id1!="NA")
personal_help_id2<-subset(personal_help_id2, personal_help_id2!="." & personal_help_id2!="NA")
personal_help_id3<-subset(personal_help_id3, personal_help_id3!="." & personal_help_id3!="NA")
personal_help_id4<-subset(personal_help_id4, personal_help_id4!="." & personal_help_id4!="NA")

# rename
names(personal_help_id1)[names(personal_help_id1) == 'covariates$personal_help_id1'] <- 'q31'
names(personal_help_id2)[names(personal_help_id2) == 'covariates$personal_help_id2'] <- 'q31'
names(personal_help_id3)[names(personal_help_id3) == 'covariates$personal_help_id3'] <- 'q31'
names(personal_help_id4)[names(personal_help_id4) == 'covariates$personal_help_id4'] <- 'q31'


# append
total <- rbind(personal_help_id1, personal_help_id2,personal_help_id3,personal_help_id4)


#aggregate
total$value <- 1
help_id<-aggregate(total$value, by=list(help=total$q31), FUN=sum)

# plot 
fig_help_id<-help_id %>%
    ggplot() +
    geom_bar(aes(x=reorder(help, x),y = (x)/sum(x)),stat="identity") +
    scale_y_continuous(labels=scales::percent) +
    coord_flip() + 
    labs(y='Number',x='Actor')+
    theme_bw()+
    ggtitle("From whom did you receive help?")

MP_satisfied <- covariates %>%
  dplyr::filter(satisfied_personal_help != ".") %>%
  dplyr::filter(personal_help_id1 == "MP",
                personal_help_id2 == ".",
                personal_help_id3 == ".",
                personal_help_id4 == ".") %>%
    dplyr::count(satisfied_personal_help) %>%
    dplyr::rename(sum_group = n) %>%
    ggplot() + 
    geom_bar(aes(x=satisfied_personal_help,y = (sum_group)/sum(sum_group)*100),stat="identity",width=0.4) +
    theme_bw()+
    labs(y='Number of respondents in %',x='How satisfied were you with the help you received?')+
    ggtitle("Subset: Respondents who recieved help only from MP")

help_id <- ggarrange(fig_help_id,MP_satisfied,
                    ncol = 2, nrow = 1)
help_id

ggsave("2_output/fig_A6_help_id.pdf", width = 12, height = 4, units = "in")

```

## Table A1: Summary Statistics Conjoint Experiment

```{r summary_stats_table1, results="asis", message=FALSE, echo=FALSE,warning=FALSE}

######################################
# summary Statistics
######################################
df_long <- readRDS("1_input/conjoint.rds")
st(df_long)
st(df_long, out = 'latex', file = "2_output/tab_A1.tex")

```

## Table A2: Summary Statistics
```{r summary_stats_table2, results="asis", message=FALSE, echo=FALSE,warning=FALSE}
covariates <- readRDS("1_input/covariates.rds") %>%
  select(-starts_with("T_"))

st(covariates)
st(covariates, out = 'latex', file = "2_output/tab_A2.tex")


```

## Table A4: Main Results
```{r , results="asis", message=FALSE, echo=FALSE,warning=FALSE}

main<-lm_robust(Choice ~Effort_prevention+Effort_relief+Effective_prevention+Effective_relief+Ask+EQ+Honesty,
                  data=df_long,clusters = CaseID, se_type = "stata")
# Update term names in the model object
main$term <- recode(
  main$term,
  `(Intercept)` = "Intercept",
  `Effort_preventionPrevention Effort` = "Preparedness Effort",
  `Effort_reliefRelief Effort` = "Relief Effort",
  `Effective_preventionPrevention Effective` = "Preparedness Effective",
  `Effective_reliefRelief Effective` = "Relief Effective",
  `AskAsk Relief` = "Relief Ask",
  `EQVisits` = "Visits",
  `HonestyCorruption` = "Corruption",
  `HonestyVote Buying` = "Vote Buying"
)
msummary(main, 
         title = 'Main Results', 
         output = "2_output/tab_A4_main_results.tex", 
         stars = TRUE, 
         fmt = 2)

msummary(main, title = 'Main Results', output = "html", stars = TRUE, fmt = 2)

```

## Figure A9: Marginal Means
```{r , results="asis", message=FALSE, echo=FALSE}
means<-mm(df_long, Choice ~Effort_prevention+Effort_relief+Effective_prevention+Effective_relief+Ask+EQ+Honesty ,id = ~ CaseID, h0 = 0.5)


means$group[means$level=="No Prevention Effort"] <- "Effort"
means$group[means$level=="Prevention Effort"] <- "Effort"
means$group[means$level=="No Relief Effort"] <- "Effort"
means$group[means$level=="Relief Effort"] <- "Effort"

means$group[means$level=="Prevention Not Effective"] <- "Effevtive"
means$group[means$level=="Prevention Effective"] <- "Effevtive"
means$group[means$level=="No Relief"] <- "Effevtive"
means$group[means$level=="Relief Effective"] <- "Effevtive"

means$group[means$level=="Not Ask Relief"] <- "Other"
means$group[means$level=="Ask Relief"] <- "Other"
means$group[means$level=="No Visits"] <- "Other"
means$group[means$level=="Visits"] <- "Other"

means$group[means$level=="No Corruption"] <- "Other"
means$group[means$level=="Corruption"] <- "Other"
means$group[means$level=="Vote Buying"] <- "Other"

fig_means<-means%>% 
  mutate(level = fct_reorder(level, desc(-estimate))) %>%
  ggplot( aes(x = estimate, y = level, color = as.factor(group))) +
  geom_point(size = 2.5) + 
  geom_errorbarh(aes(y = level, 
                     xmin = estimate - 1.96*std.error, 
                     xmax = estimate + 1.96*std.error),
                 size=0.5, alpha = 0.5, height = 0.2) +
  scale_colour_viridis_d(name = 'group') + 
  theme_bw() +
  ggtitle("Marginal Means")+
  theme(legend.position = "bottom")+
  theme(axis.title.y=element_blank()) +
  theme(axis.text=element_text(size = 12)) +
  theme(axis.text.x = element_text(size = 12))+
  geom_vline(xintercept = 0.5, linetype="dotted", 
                color = "red")

fig_means


ggsave("2_output/fig_A9_mm.pdf", width = 6, height = 4, units = "in")
```


## Table A5: Main Results, Linear Hypothesis 1

```{r , results="asis", message=FALSE, echo=FALSE}
lhro1 <- lh_robust(Choice ~Effort_prevention+Effort_relief+Effective_prevention+Effective_relief+Ask+EQ+Honesty, data = df_long, linear_hypothesis = "Effort_preventionPrevention Effort = Effort_reliefRelief Effort")

# Rename terms in the lm_robust part
lhro1$lm_robust$term <- recode(
  lhro1$lm_robust$term,
  `(Intercept)` = "Intercept",
  `Effort_preventionPrevention Effort` = "Preparedness Effort",
  `Effort_reliefRelief Effort` = "Relief Effort",
  `Effective_preventionPrevention Effective` = "Preparedness Effective",
  `Effective_reliefRelief Effective` = "Relief Effective",
  `AskAsk Relief` = "Relief Ask",
  `EQVisits` = "Visits",
  `HonestyCorruption` = "Corruption",
  `HonestyVote Buying` = "Vote Buying"
)

# Rename terms in the linear hypothesis (lh) part
lhro1$lh$term <- recode(
  lhro1$lh$term,
  `Effort_preventionPrevention Effort = Effort_reliefRelief Effort` = "Preparedness Effort = Relief Effort"
)

msummary(lhro1, 
         title = 'Main Results, Linear Hypothesis 1', 
         output = "2_output/tab_A5_linear_hypothesis_1.tex", 
         stars = TRUE, 
         fmt = 2)

msummary(lhro1, title = 'Main Results, Linear hypothesis 1',output = "html",stars = TRUE,fmt = 2)


```

## Table A6: Main Results, Linear Hypothesis 2
```{r , results="asis", message=FALSE, echo=FALSE,warning=FALSE}
lhro2 <- lh_robust(Choice ~Effort_prevention+Effort_relief+Effective_prevention+Effective_relief+Ask+EQ+Honesty, data = df_long, linear_hypothesis = "Effective_preventionPrevention Effective  = Effective_reliefRelief Effective ")


# Rename terms in the lm_robust part
lhro2$lm_robust$term <- recode(
  lhro2$lm_robust$term,
  `(Intercept)` = "Intercept",
  `Effort_preventionPrevention Effort` = "Preparedness Effort",
  `Effort_reliefRelief Effort` = "Relief Effort",
  `Effective_preventionPrevention Effective` = "Preparedness Effective",
  `Effective_reliefRelief Effective` = "Relief Effective",
  `AskAsk Relief` = "Relief Ask",
  `EQVisits` = "Visits",
  `HonestyCorruption` = "Corruption",
  `HonestyVote Buying` = "Vote Buying"
)

# Rename terms in the linear hypothesis (lh) part
lhro2$lh$term <- recode(
  lhro2$lh$term,
  `Effective_preventionPrevention Effective  = Effective_reliefRelief Effective ` = "Preparedness Effective = Relief Effective"
)

msummary(lhro2, 
         title = 'Main Results, Linear Hypothesis 2', 
         output = "2_output/tab_A6_linear_hypothesis_2.tex", 
         stars = TRUE, 
         fmt = 2)

msummary(lhro2, title = 'Main Results, Linear hypothesis 2',output = "html",stars = TRUE,fmt = 2)


```

## Table A7: Main Results with Interactions

```{r , results="asis", message=FALSE, echo=FALSE,warning=FALSE}
interaction_model <- lm_robust(
  Choice ~ Effort_prevention+Effort_relief +
    Effective_prevention + Effective_relief +
    Honesty+
    Effective_prevention* Honesty+
    Effort_relief*Effective_relief+
    Effective_prevention* Effective_relief+
    Effort_prevention * Effort_relief+
    Effort_prevention * Effective_prevention+
    Effort_prevention * Honesty +
    Effective_relief*Honesty,
    data = df_long,
  clusters = CaseID,
  se_type = "stata"
)

# Rename terms in the interaction model
# Rename terms in the interaction model
interaction_model$term <- recode(
  interaction_model$term,
  `(Intercept)` = "Intercept",
  `Effort_preventionPrevention Effort` = "Preparedness Effort",
  `Effort_reliefRelief Effort` = "Relief Effort",
  `Effective_preventionPrevention Effective` = "Preparedness Effective",
  `Effective_reliefRelief Effective` = "Relief Effective",
  `HonestyCorruption` = "Corruption",
  `HonestyVote Buying` = "Vote Buying",
  `Effort_reliefRelief Effort:Effective_reliefRelief Effective` = "Relief Effort × Relief Effective",
  `Effective_preventionPrevention Effective:Effective_reliefRelief Effective` = "Preparedness Effective × Relief Effective",
  `Effort_preventionPrevention Effort:Effort_reliefRelief Effort` = "Preparedness Effort × Relief Effort",
  `Effort_preventionPrevention Effort:Effective_preventionPrevention Effective` = "Preparedness Effort × Preparedness Effective",
  `Effort_preventionPrevention Effort:HonestyCorruption` = "Preparedness Effort × Corruption",
  `Effort_preventionPrevention Effort:HonestyVote Buying` = "Preparedness Effort × Vote Buying",
  `Effective_preventionPrevention Effective:HonestyCorruption` = "Preparedness Effective × Corruption",
  `Effective_preventionPrevention Effective:HonestyVote Buying` = "Preparedness Effective × Vote Buying",
  `Effective_reliefRelief Effective:HonestyCorruption` = "Relief Effective × Corruption",
  `Effective_reliefRelief Effective:HonestyVote Buying` = "Relief Effective × Vote Buying"
)

# Display the table in QMD and save as LaTeX
msummary(
  interaction_model,
  title = 'Table A7: Main Results with Interactions', 
  output = "2_output/tab_A7_interactions.tex", 
  stars = TRUE,
  fmt = 2
)
msummary(interaction_model, title = 'Table A7: Main Results with Interactions',output = "html",stars = TRUE,fmt = 2)
         
```



## Figure A10: Heterogeneity of AMCE by Number of Contest
```{r , results="asis", message=FALSE, echo=FALSE,warning=FALSE}
  
  df_long_main<-df_long %>% dplyr::rename(
                          EffortPrevention_ = Effort_prevention,
                          EffortRelief_ = Effort_relief,
                          EffectivePrevention_ = Effective_prevention,
                          EffectiveRelief_ = Effective_relief,
                          Ask_ = Ask,
                          EQ_ = EQ,
                          Honesty_ = Honesty
                          )
  # Create an empty data frame to store the results
  final_results <- data.frame()

  # Get unique contests
  unique_contests <- unique(df_long_main$contest)

  # Loop through each unique contest
  for (contest_value in unique_contests) {
  # Create a subset of the dataframe for the current contest value
  contest_data <- df_long_main %>% filter(contest == contest_value)
  
  # Perform the analysis for the current contest
  res <- lm_robust(Choice ~EffortPrevention_+EffortRelief_+EffectiveRelief_+EffectivePrevention_+Ask_+EQ_+Honesty_,
                   data = contest_data, clusters = CaseID, se_type = "stata") %>%
    tidy_coefs_with_ref() %>%
    dplyr::filter(!(variable == '(Intercept)_')) %>%
    dplyr::mutate(level = paste(level, "   ")) %>%
    dplyr::mutate(estimate = tidyr::replace_na(estimate, 0)) %>%
    dplyr::mutate(
      variable = dplyr::recode(variable, 
                               "Ask_"="E. Ask_", 
                               "Honesty_" = "G. Honesty_",
                               "EQ_" = "F. EQ_",
                               "EffectivePrevention_" = "C. EffectivePrevention_",
                               "EffectiveRelief_" = "D. EffectiveRelief_",
                               "EffortPrevention_" = "A. EffortPrevention_",
                               "EffortRelief_" = "B. EffortRelief_"
      )) %>%
    # add empty raw
    as_tibble() %>%
    group_by(variable) %>% 
    group_modify(~ add_row(.x,.before=0)) %>%
    # reorder
    dplyr::mutate(level = fct_reorder(level, estimate,.na_rm = TRUE))  %>% 
    dplyr::mutate(level = dplyr::coalesce(level,variable))   %>% 
    dplyr::mutate(
      level = dplyr::recode(level, 
                            "A. EffortPrevention_" = "Effort Preparedness",
                            "B. EffortRelief_" = "Effort Relief",
                            "C. EffectivePrevention_" = "Effective Preparedness",
                            "D. EffectiveRelief_" = "Effective Relief",
                            "E. Ask_"="Ask Help", 
                            "F. EQ_" = "Visits",
                            "G. Honesty_" = "Corruption"
      )) %>%
    add_column(indicator = 0) %>%
    dplyr::mutate(
      indicator = replace(indicator, variable == "A. EffortPrevention_", "Effort"),
      indicator = replace(indicator, variable == "B. EffortRelief_", "Effort"),
      indicator = replace(indicator, variable == "C. EffectivePrevention_", "Effective"),
      indicator = replace(indicator, variable == "D. EffectiveRelief_", "Effective"),
      indicator = replace(indicator, variable == "E. Ask_", "Other"),
      indicator = replace(indicator, variable == "F. EQ_", "Other"),
      indicator = replace(indicator, variable == "G. Honesty_", "Other")
    ) %>%
    dplyr::mutate(
      level = dplyr::recode(level, 
                            "Prevention Effective    " = "Yes    ", 
                            "Prevention Not Effective    " = "No    ",
                            "Relief Effective    " = " Yes    ",
                            "No Relief    " = " No    ",
                            "Prevention Effort    " = "High    ",
                            "No Prevention Effort    " = "Low    ",
                            "Relief Effort    " = " High   ",
                            "No Relief Effort    " = " Low    ",
                            "Ask Relief    " = "  Yes    ",
                            "Not Ask Relief    " = "  No    ",
                            "No Visits    " = "   No    ",
                            "Visits    " = "   Yes    "
      )
    ) %>%
    mutate(level = factor(level, levels = level)) %>%  
    mutate(level = fct_rev(as_factor(level)))
    # Add a column indicating the contest for each observation
    res$contest <- paste("Contest", contest_value)
  
  # Bind the results for the current contest to the final_results data frame
  final_results <- bind_rows(final_results, res)
  }


  bold.ticks <- c("Ask Help","Corruption","Visits","Effective Preparedness","Effective Relief","Effort Preparedness","Effort Relief")
  bold.labels <- ifelse(levels(final_results$level) %in% bold.ticks, yes = "bold", no = "plain")

  fig_contest<-final_results %>%  
  ggplot( aes(x = estimate, y =  level,color = indicator)) +
    geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = level, 
                     xmin = conf.low, 
                     xmax = conf.high),
                 size=0.5, alpha = 1, height = 0.2) +
  scale_colour_colorblind()+
  facet_wrap(vars(as.factor(contest)))+
  theme_bw() +
  ggtitle("AMCEs")+
  xlab("Effect on Probability of Support for Candidate") +
  ylab("Disaster Policy Choice")+
  theme(axis.text=element_text(size = 12)) +
  theme(strip.text.x = element_text(size =12))+
  theme(axis.text.y = element_text(face=bold.labels))+
  theme(axis.title = element_text(size=12,face="bold",vjust=-2.5))+
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")+
  theme(legend.position = "right")

  fig_contest
  
  ggsave("2_output/fig_A10_contest.pdf", width = 14, height = 8, units = "in")
  
  
```

## Figure A11: Distance to Flood (2015) in 2018.
```{r flood_distance, message=FALSE, warning=FALSE, echo=FALSE,out.width="60%",fig.align="center",fig.cap=c("Distance of survey respindents to Flood (2015) in 2018 (in meters)")}

  df_long_main <- left_join(df_long,covariates, by = c("CaseID")) 

  df_long_main<-df_long_main %>% dplyr::rename(
                          EffortPrevention_ = Effort_prevention,
                          EffortRelief_ = Effort_relief,
                          EffectivePrevention_ = Effective_prevention,
                          EffectiveRelief_ = Effective_relief,
                          Ask_ = Ask,
                          EQ_ = EQ,
                          Honesty_ = Honesty
                          )


  df_distance <- df_long_main  %>% 
  dplyr::filter(life2015 == "2")%>% 
  dplyr::mutate(distance_cat = case_when(
                              distance_flood >=0 & distance_flood <=1000 ~ "A. 0-1000m",
                              distance_flood >=1000 & distance_flood <=2500 ~ "B. 1000-2500m",
                              distance_flood >=2500 & distance_flood <=5000 ~ "C. 2500-5000m",
                              distance_flood >=5000 & distance_flood <=10000 ~ "D. 5000-10000m",
                              distance_flood >=10000 & distance_flood <=15000 ~ "E. 10000-15000m",
                              distance_flood >=15000 ~ "F. 15000-20000m"))
  # plot distance
  flood_distance2018<-ggplot() +
  geom_histogram(df_distance, mapping=aes(as.numeric(distance_flood)),alpha = 0.7)+
  labs(y='Density',x='Distance to Flood (2015) in 2018 (meters)')+
  theme_bw() 
  flood_distance2018

  ggsave("2_output/fig_A11_flood_distance2018.pdf", width = 6, height = 3, units = "in")

```
## Figure 2: AMCE’s by by distance to 2015 flood

```{r , results="asis", message=FALSE, echo=FALSE}
  
  # Create an empty data frame to store the results
  final_results <- data.frame()

  # Get unique contests
  unique_distance <- unique(df_distance$distance_cat)

  # Loop through each unique contest
  for (distance_value in unique_distance) {
  # Create a subset of the dataframe for the current distance_cat value
  distance_data <- df_distance %>% filter(distance_cat == distance_value)
  
  # Perform the analysis for the current distance_cat
  res <- lm_robust(Choice ~EffortPrevention_+EffortRelief_+EffectiveRelief_+EffectivePrevention_+Ask_+EQ_+Honesty_,
                   data = distance_data, clusters = CaseID, se_type = "stata") %>%
    tidy_coefs_with_ref() %>%
    dplyr::filter(!(variable == '(Intercept)_')) %>%
    dplyr::mutate(level = paste(level, "   ")) %>%
    dplyr::mutate(estimate = tidyr::replace_na(estimate, 0)) %>%
    dplyr::mutate(
      variable = dplyr::recode(variable, 
                               "Ask_"="E. Ask_", 
                               "Honesty_" = "G. Honesty_",
                               "EQ_" = "F. EQ_",
                               "EffectivePrevention_" = "C. EffectivePrevention_",
                               "EffectiveRelief_" = "D. EffectiveRelief_",
                               "EffortPrevention_" = "A. EffortPrevention_",
                               "EffortRelief_" = "B. EffortRelief_"
      )) %>%
    # add empty raw
    as_tibble() %>%
    group_by(variable) %>% 
    group_modify(~ add_row(.x,.before=0)) %>%
    # reorder
    dplyr::mutate(level = fct_reorder(level, estimate,.na_rm = TRUE))  %>% 
    dplyr::mutate(level = dplyr::coalesce(level,variable))   %>% 
    dplyr::mutate(
      level = dplyr::recode(level, 
                            "A. EffortPrevention_" = "Effort Preparedness",
                            "B. EffortRelief_" = "Effort Relief",
                            "C. EffectivePrevention_" = "Effective Preparedness",
                            "D. EffectiveRelief_" = "Effective Relief",
                            "E. Ask_"="Ask Help", 
                            "F. EQ_" = "Visits",
                            "G. Honesty_" = "Corruption"
      )) %>%
    add_column(indicator = 0) %>%
    dplyr::mutate(
      indicator = replace(indicator, variable == "A. EffortPrevention_", "Effort"),
      indicator = replace(indicator, variable == "B. EffortRelief_", "Effort"),
      indicator = replace(indicator, variable == "C. EffectivePrevention_", "Effective"),
      indicator = replace(indicator, variable == "D. EffectiveRelief_", "Effective"),
      indicator = replace(indicator, variable == "E. Ask_", "Other"),
      indicator = replace(indicator, variable == "F. EQ_", "Other"),
      indicator = replace(indicator, variable == "G. Honesty_", "Other")
    ) %>%
    dplyr::mutate(
      level = dplyr::recode(level, 
                            "Prevention Effective    " = "Yes    ", 
                            "Prevention Not Effective    " = "No    ",
                            "Relief Effective    " = " Yes    ",
                            "No Relief    " = " No    ",
                            "Prevention Effort    " = "High    ",
                            "No Prevention Effort    " = "Low    ",
                            "Relief Effort    " = " High   ",
                            "No Relief Effort    " = " Low    ",
                            "Ask Relief    " = "  Yes    ",
                            "Not Ask Relief    " = "  No    ",
                            "No Visits    " = "   No    ",
                            "Visits    " = "   Yes    "
      )
    ) %>%
    mutate(level = factor(level, levels = level)) %>%  
    mutate(level = fct_rev(as_factor(level)))
    # Add a column indicating the contest for each observation
    res$distance_cat <- paste("", distance_value)
  
  # Bind the results for the current contest to the final_results data frame
  final_results <- bind_rows(final_results, res)
  }


  bold.ticks <- c("Ask Help","Corruption","Visits","Effective Preparedness","Effective Relief","Effort Preparedness","Effort Relief")
  bold.labels <- ifelse(levels(final_results$level) %in% bold.ticks, yes = "bold", no = "plain")

  fig_distance_cat<-final_results %>%  
  ggplot( aes(x = estimate, y =  level,color = indicator)) +
    geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = level, 
                     xmin = conf.low, 
                     xmax = conf.high),
                 size=0.5, alpha = 1, height = 0.2) +
  scale_colour_colorblind()+
  facet_wrap(vars(as.factor(distance_cat)))+
  theme_bw() +
  ggtitle("AMCEs")+
  xlab("Effect on Probability of Support for Candidate") +
  ylab("Disaster Policy Choice")+
  theme(axis.text=element_text(size = 12)) +
  theme(strip.text.x = element_text(size =12))+
  theme(axis.text.y = element_text(face=bold.labels))+
  theme(axis.title = element_text(size=12,face="bold",vjust=-2.5))+
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")+
  theme(legend.position = "bottom")
  fig_distance_cat
  
  ggsave("2_output/fig_A12_distance_cat.pdf", width = 14, height = 8, units = "in")

```


### Figure A13: Scatterplot Economic Distress (X) and Distance to Flood (Y)


```{r scattter, message=FALSE, warning=FALSE, echo=FALSE,out.width="50%",fig.align="center",fig.cap=c("Scatter plot, Distance to the 2015 flood in meters (X-Axis) and seft-reported economic hardship due to the flood (Y-axis)")}
fig_harm_distance<-ggplot(covariates, aes(x=flood_econ, y=distance_flood)) + 
  geom_point(shape=18, color="blue")+
  xlab("How badly were you economically harmed by the 2015 floods?") +
  ylab("Distance to the 2015 flood (meters)")+
  geom_smooth(method=lm, se=FALSE, linetype="dashed",
             color="darkred")+
  theme_bw()
fig_harm_distance

  ggsave("2_output/fig_A13_harm_distance.pdf", width = 8, height = 5, units = "in")

```

### Figure A14: Flood exposure
```{r flood_econ, message=FALSE, warning=FALSE, echo=FALSE,out.width="70%",fig.align="center",message=FALSE, results='hide',fig.cap=c("Distribution: seft-reported economic hardship due to the flood")}
  covariates<-covariates %>%
  dplyr::mutate(
  # make binary and character
  flood_econ_bin = case_when(flood_econ <=3 ~ 0,
                             flood_econ >3 ~ 1
  ),
  flood_econ_char = case_when(
                              flood_econ ==1 ~ "Not at all",
                              flood_econ ==2 ~ "Just mildly",
                              flood_econ ==3 ~ "Somewhat",
                              flood_econ ==4 ~ "Very badly",
                              flood_econ ==5 ~ "Extremely badly")
  )
  # Order  
  manual_country_order <- c("Not at all","Just mildly","Somewhat","Very badly","Extremely badly")
  # Apply the manual order 
  covariates$flood_econ_char <- factor(covariates$flood_econ_char, levels = manual_country_order)  

  flood_econ<-covariates %>%
  ggplot(aes(flood_econ_char)) +
  geom_bar() +
  geom_text(aes(label = scales::percent(..prop..), group = 1), stat= "count")+ 
  scale_colour_viridis_d(name = 'S_preference') + 
  theme_bw() +
  xlab("How badly were you economically harmed by the 2015 floods?") +
  ylab("Number of observations")+
  theme(axis.text=element_text(size = 12)) +
  theme(strip.text.x = element_text(size =12))+
  theme(legend.position = "none")
  flood_econ
  
  ggsave("2_output/fig_A14_flood_econ.pdf", width = 8, height = 4, units = "in")

```

##  Figure A15: Difference in AMCEs and Marginal Means

```{r, message=FALSE, warning=FALSE, echo=FALSE,out.width="100%",fig.align="center",message=FALSE, results='hide'}
  df_long_main <- left_join(df_long,covariates, by = c("CaseID")) 

  df_long_main<-df_long_main %>%
  dplyr::mutate(
  # make binary and character
  flood_econ_bin = case_when(flood_econ <=3 ~ 0,
                             flood_econ >3 ~ 1
  ))
  
economic<-lm_robust(Choice ~Effort_prevention+Effort_relief+Effective_prevention+Effective_relief+Ask+EQ+Honesty+flood_econ_bin+Effort_prevention*flood_econ_bin+
                   Effort_relief*flood_econ_bin+
                   Effective_prevention*flood_econ_bin+
                   Effective_relief*flood_econ_bin+
                   Ask*flood_econ_bin+
                   EQ*flood_econ_bin+
                   Honesty*flood_econ_bin,
                   data=df_long_main,clusters = CaseID) %>% tidy() %>% 
  dplyr::filter(grepl(':', term)) %>% 
  dplyr::mutate(term=dplyr::recode(term, 
                     "Effort_preventionPrevention Effort:flood_econ_bin"="Preparedness Effort", 
                     "Effort_reliefRelief Effort:flood_econ_bin" = "Relief Effort",
                     "Effective_preventionPrevention Effective:flood_econ_bin" = "Preparedness Effevtive",
                     "Effective_reliefRelief Effective:flood_econ_bin" = "Relief Effevtive",
                     "AskAsk Relief:flood_econ_bin" = "Ask Help",
                     "EQVisits:flood_econ_bin" = "Visits",
                     "HonestyCorruption:flood_econ_bin" = "Corruption",
                     "HonestyVote Buying:flood_econ_bin" = "Vote Buying"
                     )
         )%>% 
  dplyr::mutate(
  group = case_when(
  term == "Preparedness Effort"  ~ "Effort",
  term == "Relief Effort"  ~ "Effort",
  term == "Preparedness Effevtive"  ~ "Effevtive",
  term == "Relief Effevtive"  ~ "Effevtive",
  term == "Ask Help"  ~ "Other",
  term == "Visits"  ~ "Other",
  term == "Corruption"  ~ "Other",
  term == "Vote Buying"  ~ "Other"
  )
  )

f7<-economic%>% 
  mutate(term = fct_reorder(term, desc(-estimate))) %>%
  ggplot( aes(x = estimate, y = term)) +#color = as.factor(group)
      geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = term, 
                     xmin =conf.low, 
                     xmax = conf.high),
                 size=0.5, alpha = 0.5, height = 0.2) +
  #scale_colour_viridis_d(name = 'Group') + 
  theme_bw() +
  ggtitle("difference in AMCEs (low vs. high losses)")+
  scale_x_continuous(limits = c(-0.27, 0.2))+
  theme(axis.title.y=element_blank()) +
  theme(axis.text=element_text(size = 12)) +
  theme(axis.text.x = element_text(size = 12)) +
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")+
  xlab("Effect on Probability of Support for Candidate") 
# +
    #theme(legend.position="bottom")

  df_long_main$flood_econ_bin <- as.factor(df_long_main$flood_econ_bin) 


f9<-cj(df_long_main, Choice~Effort_prevention+Effort_relief+Effective_prevention+Effective_relief+Ask+EQ+Honesty, id = ~CaseID, estimate = "mm_differences", by = ~flood_econ_bin)%>% 
  dplyr::mutate(level = fct_reorder(level, desc(-estimate))) %>%
  ggplot( aes(x = estimate, y = level)) +# color = as.factor(group)
      geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = level, 
                     xmin =lower, 
                     xmax = upper),
                 size=0.5, alpha = 0.5, height = 0.2) +
    scale_colour_viridis_d(name = 'Group') + 
  theme_bw() +
  ggtitle("difference in marginal means (low vs. high losses)")+
  scale_x_continuous(limits = c(-0.27, 0.2))+
  theme(axis.title.y=element_blank()) +
  theme(axis.text=element_text(size = 12)) +
  theme(axis.text.x = element_text(size = 12)) +
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")+
  xlab("Effect on Probability of Support for Candidate") 


fig_econ  <- ggarrange(f7,f9,
          ncol = 2, nrow = 1)

fig_econ
ggsave("2_output/fig_A15_econ.pdf", width = 12, height = 4, units = "in")

```

## Figure A16: Balance Tests

```{r, message=FALSE, warning=FALSE, echo=FALSE,out.width="100%",fig.align="center",message=FALSE, results='hide'}
###################
# Balance Check
###################
# subset to farmer (only here the prime should work)
total1 <- subset(df_long_main, farmer==1)
head(total1)
total1$disaster_prime[total1$disaster_prime=="control"] <- 0
total1$disaster_prime[total1$disaster_prime=="treatment"] <- 1
#total1$gender[total1$gender=="Female"] <- 0
#total1$gender[total1$gender=="Male"] <- 1
total1<-total1[!(total1$gender>1),]
summary(total1$gender)


# select variables 
myvars <- c("disaster_prime", "education", "income", "age","gender","interested_politics","trust_MP")
newdata <- total1[myvars]


# Modify each plot to include a legend with a common title
education <- ggplot() + 
  geom_density(data = newdata, aes(x = education, fill = factor(disaster_prime, labels = c("Control", "Treatment"))), 
               alpha = 0.2, adjust = 2) + 
  xlab("Education") +
  ylab("Density") +
  theme_classic() + 
  scale_fill_grey() +
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 14)) +
  labs(fill = "Treatment Group")

age <- ggplot() + 
  geom_density(data = newdata, aes(x = age, fill = factor(disaster_prime, labels = c("Control", "Treatment"))), 
               alpha = 0.2, adjust = 2) + 
  xlab("Age") +
  ylab("Density") +
  theme_classic() + 
  scale_fill_grey() +
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 14)) +
  labs(fill = "Treatment Group")

income <- ggplot() + 
  geom_density(data = newdata, aes(x = income, fill = factor(disaster_prime, labels = c("Control", "Treatment"))), 
               alpha = 0.2, adjust = 2) + 
  xlab("Income") +
  ylab("Density") +
  theme_classic() + 
  scale_fill_grey() +
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 14)) +
  labs(fill = "Treatment Group")

gender <- ggplot() + 
  geom_density(data = newdata, aes(x = gender, fill = factor(disaster_prime, labels = c("Control", "Treatment"))), 
               alpha = 0.2, adjust = 2) + 
  xlab("Gender") +
  ylab("Density") +
  theme_classic() + 
  scale_fill_grey() +
  theme(axis.text = element_text(size = 12),
        axis.title = element_text(size = 14)) +
  labs(fill = "Treatment Group")

# Combine the plots into a grid with a common legend
balance <- ggarrange(education, income, age, gender,
                     ncol = 2, nrow = 2, common.legend = TRUE, legend = "bottom")

# Display the combined plot
balance

ggsave("2_output/fig_A16_balance.pdf", width = 8, height = 4, units = "in")

```

## Table A8: Manipulation Check

```{r ,results = 'asis'}
covariates <- readRDS("1_input/covariates.rds")

manipulation<-lm_robust(manipulation ~disaster_prime,data=covariates)
resp_MP<-lm_robust(T_q37_4 ~disaster_prime,data=covariates)
resp_MP2<-lm_robust(T_q38_4 ~disaster_prime,data=covariates)
worried<-lm_robust(worried ~disaster_prime,data=covariates)
hopeless<-lm_robust(T_q24_4 ~disaster_prime,data=covariates)

# Table 
tex_table <- texreg(
  list(manipulation, resp_MP, resp_MP2, worried, hopeless),
  include.ci = FALSE,
  stars = c(0.001, 0.01, 0.05, 0.1),
  caption = "Manipulation Check",
  custom.model.names = c("Financial Worry", "MP Responsible Relief", 
                         "MP Responsible Preparedness", "Flood Worry", "Hopeless")
)

# Write the LaTeX table to a file
write(tex_table, file = "2_output/tab_A8_manipulation_check.tex")

htmlreg(
  list(manipulation, resp_MP, resp_MP2, worried, hopeless),
  include.ci = FALSE,
  stars = c(0.001, 0.01, 0.05, 0.1),
  caption = "Manipulation Check",
  custom.model.names = c("Financial Worry", "MP Responsible Relief", 
                         "MP Responsible Preparedness", "Flood Worry", "Hopeless")
)


```


## Figure A17: Difference in AMCEs of Attributes on Candidate Support, by disaster prime

```{r ,results = 'asis',warning=FALSE}
  df_long_main <- left_join(df_long,covariates, by = c("CaseID")) 

# Prime
prime<-lm_robust(Choice ~Effort_prevention+Effort_relief+Effective_prevention+Effective_relief+Ask+EQ+Honesty+disaster_prime+Effort_prevention*disaster_prime+
                   Effort_relief*disaster_prime+
                   Effective_prevention*disaster_prime+
                   Effective_relief*disaster_prime+
                   Ask*disaster_prime+
                   EQ*disaster_prime+
                   Honesty*disaster_prime,
                   data=df_long_main,clusters = CaseID) %>% tidy() %>% 
  dplyr::filter(grepl(':', term)) %>% 
  dplyr::mutate(term=dplyr::recode(term, 
                     "Effort_preventionPrevention Effort:disaster_primetreatment"="Preparedness Effort", 
                     "Effort_reliefRelief Effort:disaster_primetreatment" = "Relief Effort",
                     "Effective_preventionPrevention Effective:disaster_primetreatment" = "Preparedness Effevtive",
                     "Effective_reliefRelief Effective:disaster_primetreatment" = "Relief Effevtive",
                     "AskAsk Relief:disaster_primetreatment" = "Ask Help",
                     "EQVisits:disaster_primetreatment" = "Visits",
                     "HonestyCorruption:disaster_primetreatment" = "Corruption",
                     "HonestyVote Buying:disaster_primetreatment" = "Vote Buying"
                     )
         )%>% 
  dplyr::mutate(
  group = case_when(
  term == "Preparedness Effort"  ~ "Effort",
  term == "Relief Effort"  ~ "Effort",
  term == "Preparedness Effevtive"  ~ "Effevtive",
  term == "Relief Effevtive"  ~ "Effevtive",
  term == "Ask Help"  ~ "Other",
  term == "Visits"  ~ "Other",
  term == "Corruption"  ~ "Other",
  term == "Vote Buying"  ~ "Other"
  )
  )


fig_prime<-prime%>% 
  mutate(term = fct_reorder(term, desc(-estimate))) %>%
  ggplot( aes(x = estimate, y = term)) +# color = as.factor(group)
      geom_point(size = 2.5) + 
    geom_errorbarh(aes(y = term, 
                     xmin =conf.low, 
                     xmax = conf.high),
                 size=0.5, alpha = 0.5, height = 0.2) +
  #scale_colour_viridis_d(name = 'Group') + 
  theme_bw() +
  ggtitle("Disaster Prime")+
  scale_x_continuous(limits = c(-0.27, 0.2))+
  theme(axis.title.y=element_blank()) +
  theme(axis.text=element_text(size = 12)) +
  theme(axis.text.x = element_text(size = 12)) +
  geom_vline(xintercept = 0, linetype="dotted", 
                color = "red")
fig_prime
ggsave("2_output/fig_A17_prime.pdf", width = 8, height = 4, units = "in")


```


## Table A9: Effect of Economic Losses, Full Models, with Controls

```{r ,results = 'asis',warning=FALSE}

  df_long_main <- left_join(df_long,covariates, by = c("CaseID")) 

  df_long_main<-df_long_main %>%
  dplyr::mutate(
  # make binary and character
  flood_econ_bin = case_when(flood_econ <=3 ~ 0,
                             flood_econ >3 ~ 1
  ))
  
# Economic Distress
econ_distress1<-lm_robust(Choice ~Effort_prevention+Effort_relief+Ask+Effective_prevention+Effective_relief+EQ+Honesty+flood_econ_bin+Effort_prevention*flood_econ_bin+Effort_relief*flood_econ_bin+Ask*flood_econ_bin+Effective_prevention*flood_econ_bin+Effective_relief*flood_econ_bin+EQ*flood_econ_bin+Honesty*flood_econ_bin,data=df_long_main,clusters = CaseID)
# controlling for poverty
econ_distress2<-lm_robust(Choice ~Effort_prevention+Effort_relief+Ask+Effective_prevention+Effective_relief+EQ+Honesty+flood_econ_bin+Effort_prevention*flood_econ_bin+Effort_relief*flood_econ_bin+Ask*flood_econ_bin+Effective_prevention*flood_econ_bin+Effective_relief*flood_econ_bin+EQ*flood_econ_bin+Honesty*flood_econ_bin+poverty,data=df_long_main,clusters = CaseID)
# controlling for help recieved
econ_distress3<-lm_robust(Choice ~Effort_prevention+Effort_relief+Ask+Effective_prevention+Effective_relief+EQ+Honesty+flood_econ_bin+Effort_prevention*flood_econ_bin+Effort_relief*flood_econ_bin+Ask*flood_econ_bin+Effective_prevention*flood_econ_bin+Effective_relief*flood_econ_bin+EQ*flood_econ_bin+Honesty*flood_econ_bin+poverty,data=df_long_main,clusters = CaseID)
# controlling for satisfied with help 
econ_distress3<-lm_robust(Choice ~Effort_prevention+Effort_relief+Ask+Effective_prevention+Effective_relief+EQ+Honesty+flood_econ_bin+Effort_prevention*flood_econ_bin+Effort_relief*flood_econ_bin+Ask*flood_econ_bin+Effective_prevention*flood_econ_bin+Effective_relief*flood_econ_bin+EQ*flood_econ_bin+Honesty*flood_econ_bin+poverty,data=df_long_main,clusters = CaseID)
# standard covariates
econ_distress4<-lm_robust(Choice ~Effort_prevention+Effort_relief+Ask+Effective_prevention+Effective_relief+EQ+Honesty+flood_econ_bin+Effort_prevention*flood_econ_bin+Effort_relief*flood_econ_bin+Ask*flood_econ_bin+Effective_prevention*flood_econ_bin+Effective_relief*flood_econ_bin+EQ*flood_econ_bin+Honesty*flood_econ_bin+poverty+education+age+gender+interested_politics+trust_MP,data=df_long_main,clusters = CaseID)

econ_distress5<-lm_robust(Choice ~Effort_prevention+Effort_relief+Ask+Effective_prevention+Effective_relief+EQ+Honesty+flood_econ_bin+Effort_prevention*flood_econ_bin+Effort_relief*flood_econ_bin+Ask*flood_econ_bin+Effective_prevention*flood_econ_bin+Effective_relief*flood_econ_bin+EQ*flood_econ_bin+Honesty*flood_econ_bin+poverty+education+age+gender+interested_politics+trust_MP+worried+personal_help+poverty*flood_econ_bin,data=df_long_main,clusters = CaseID)


# Define the custom coefficient names for renaming and reordering
coef_names <- list(
  "Effort_preventionPrevention Effort:flood_econ_bin" = "Prevention Effort x Economic Losses",
  "Effort_reliefRelief Effort:flood_econ_bin" = "Relief Effort x Economic Losses",
  "AskAsk Relief:flood_econ_bin" = "Relief Ask x Economic Losses",
  "Effective_preventionPrevention Effective:flood_econ_bin" = "Prevention Effective x Economic Losses",
  "Effective_reliefRelief Effective:flood_econ_bin" = "Relief Effective x Economic Losses",
  "EQVisits:flood_econ_bin" = "Visits x Economic Losses",
  "HonestyCorruption:flood_econ_bin" = "Corruption x Economic Losses",
  "HonestyVote Buying:flood_econ_bin" = "Vote Buying x Economic Losses",
  "poverty" = "Poverty",
  "education" = "Education",
  "age" = "Age",
  "gender" = "Gender",
  "interested_politics" = "Interested Politics",
  "trust_MP" = "Trust MP",
  "worried" = "Flood Worried",
  "personal_help" = "Received Help",
   "flood_econ_bin:poverty" = "Economic Losses x Poverty"
)

# Create the regression table with renamed and selected coefficients
tab_econ_losses<-texreg(
  list(econ_distress1, econ_distress2, econ_distress3, econ_distress4, econ_distress5),
  caption = "Effect of Economic Losses, Full Models, with Controls",
  custom.coef.map = coef_names,
  include.ci = FALSE,
  stars = c(0.01, 0.05, 0.1),
  dcolumn = FALSE,
  include.adjrs = FALSE,
  include.bic = TRUE,
  include.rs = TRUE,
  include.rmse = FALSE
)

# Write the LaTeX table to a file
write(tab_econ_losses, file = "2_output/tab_A9_econ_losses.tex")

htmlreg(
  list(econ_distress1, econ_distress2, econ_distress3, econ_distress4, econ_distress5),
  caption = "Effect of Economic Losses, Full Models, with Controls",
  custom.coef.map = coef_names,
  include.ci = FALSE,
  stars = c(0.01, 0.05, 0.1),
  dcolumn = FALSE,
  include.adjrs = FALSE,
  include.bic = TRUE,
  include.rs = TRUE,
  include.rmse = FALSE
)


```

## Figure A18: Media coverage. 2017-2024, Source: GDELT
```{r,warning=FALSE,message=FALSE}
# script base on: https://rpubs.com/drkblake/1007662
library(readr)

startdate <- "20170101"
enddate <- "20231231"
########################
########################
########################
# 1. government 
########################
########################
########################


### preparedness
query <- "'government preparedness'"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_preparedness <- Volume
Volume_preparedness$preparednessVolume <- Volume_preparedness$Value


###  relief
query <- "'government relief'"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_relief <- Volume
Volume_relief$reliefVolume <- Volume_relief$Value


### Merging
Volume_total <- merge(Volume_preparedness, Volume_relief, by = "Date") 


ww_overview<-Volume_total %>%
  dplyr::select(Date,preparednessVolume,reliefVolume)%>%
  reshape2::melt( id.var='Date') %>%
  ggplot() + 
  geom_line( aes(x = Date, y = value, col=variable)) +
  scale_x_date( date_labels = "%Y",date_breaks = 'year') +
  xlab('Date') +
  ylab("coverage as % of total article count")+ 
  theme_bw()+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=12))+
  scale_color_discrete(labels = c("government preparedness", "government relief"), name = "Article phrase")+
  ggtitle("news coverage, world-wide")


############################
############################
# Malawi MI
############################
############################

###  preparedness
query <- "'government preparedness' SourceCountry:MI"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=TimelineVol&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_preparedness <- Volume
Volume_preparedness$preparednessVolume <- Volume_preparedness$Value


###  relief
query <- "'government relief' SourceCountry:MI"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=TimelineVol&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_relief <- Volume
Volume_relief$reliefVolume <- Volume_relief$Value



Volume_test <- merge(Volume_relief, Volume_preparedness, by = "Date") 


mi_overview<-Volume_test %>%
  dplyr::select(Date,preparednessVolume,reliefVolume)%>%
  reshape2::melt( id.var='Date') %>%
  ggplot() + 
  geom_line( aes(x = Date, y = value, col=variable,alpha=0.3)) +
  scale_x_date( date_labels = "%Y",date_breaks = 'year') +
  xlab('Date') +
  ylab("coverage as % of total article count")+ 
  theme_bw()+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=12))+
  scale_color_discrete(labels = c("government preparedness", "government relief"), name = "Article phrase")+
  ggtitle("news coverage, Malawi")



news_pre_rel1<- ggarrange(ww_overview,mi_overview,
                      ncol = 2, nrow = 1,common.legend = TRUE, legend="bottom")


########################
########################
########################
# 2. disaster 
########################
########################
########################
########################

### preparedness
query <- "'disaster preparedness'"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_preparedness <- Volume
Volume_preparedness$preparednessVolume <- Volume_preparedness$Value


###  relief
query <- "'disaster relief'"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_relief <- Volume
Volume_relief$reliefVolume <- Volume_relief$Value



### Merging
Volume_total <- merge(Volume_preparedness, Volume_relief, by = "Date") 


ww_overview2<-Volume_total %>%
  dplyr::select(Date,preparednessVolume,reliefVolume)%>%
  reshape2::melt( id.var='Date') %>%
  ggplot() + 
  geom_line( aes(x = Date, y = value, col=variable)) +
  scale_x_date( date_labels = "%Y",date_breaks = 'year') +
  xlab('Date') +
  ylab("coverage as % of total article count")+ 
  theme_bw()+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=12))+
  scale_color_discrete(labels = c("disaster preparedness", "disaster relief"), name = "Article phrase")+
  ggtitle("news coverage, world-wide")


############################
############################
# Malawi MI
############################
############################

###  preparedness
query <- "'disaster preparedness' SourceCountry:MI"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=TimelineVol&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_preparedness <- Volume
Volume_preparedness$preparednessVolume <- Volume_preparedness$Value


###  relief
query <- "'disaster relief' SourceCountry:MI"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=TimelineVol&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_relief <- Volume
Volume_relief$reliefVolume <- Volume_relief$Value


### natural disaster 
query <- "'natural disaster' SourceCountry:MI"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_disaster <- Volume
Volume_disaster$disasterVolume <- Volume_disaster$Value

### Merging
Volume_total <- merge(Volume_preparedness, Volume_relief, by = "Date") 



Volume_test <- merge(Volume_total, Volume_disaster, by = "Date") 


mi_overview2<-Volume_test %>%
  dplyr::select(Date,preparednessVolume,reliefVolume)%>%
  reshape2::melt( id.var='Date') %>%
  ggplot() + 
  geom_line( aes(x = Date, y = value, col=variable,alpha=0.3)) +
  scale_x_date( date_labels = "%Y",date_breaks = 'year') +
  xlab('Date') +
  ylab("coverage as % of total article count")+ 
  theme_bw()+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=12))+
  scale_color_discrete(labels = c("disaster preparedness", "disaster relief"), name = "Article phrase")+
  ggtitle("news coverage, Malawi")



news_pre_rel2<- ggarrange(ww_overview2,mi_overview2,
                         ncol = 2, nrow = 1,common.legend = TRUE, legend="bottom")



########################
########################
########################
# 3. Benchmark 
########################
########################
########################

### taxation
query <- "'taxation'"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_taxation <- Volume
Volume_taxation$taxationVolume <- Volume_taxation$Value


###  unemployment
query <- "'unemployment'"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_unemployment <- Volume
Volume_unemployment$unemploymentVolume <- Volume_unemployment$Value



### Merging
Volume_total <- merge(Volume_taxation, Volume_unemployment, by = "Date") 


ww_overview<-Volume_total %>%
  dplyr::select(Date,taxationVolume,unemploymentVolume)%>%
  reshape2::melt( id.var='Date') %>%
  ggplot() + 
  geom_line( aes(x = Date, y = value, col=variable)) +
  scale_x_date( date_labels = "%Y",date_breaks = 'year') +
  xlab('Date') +
  ylab("coverage as % of total article count")+ 
  theme_bw()+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=12))+
  scale_color_discrete(labels = c("taxation", "unemployment"), name = "Article phrase")+
  ggtitle("news coverage, world-wide")


############################
############################
# Malawi MI
############################
############################

### taxation
query <- "'taxation' SourceCountry:MI"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_taxation <- Volume
Volume_taxation$taxationVolume <- Volume_taxation$Value


###  unemployment
query <- "'unemployment' SourceCountry:MI"
#Building the Volume dataframe
vp1 <- "https://api.gdeltproject.org/api/v2/doc/doc?query="
vp2 <- "&mode=timelinevolinfo&startdatetime="
vp3 <- "000000&enddatetime="
vp4 <- "000000&format=CSV"
text_v_url <- paste0(vp1, query, vp2, startdate, vp3, enddate, vp4)
v_url <- URLencode(text_v_url)
v_url

Volume <- read_csv(v_url)

Volume$Date <- as.Date(Volume$Date, "%Y-%m-%d")
Volume_unemployment <- Volume
Volume_unemployment$unemploymentVolume <- Volume_unemployment$Value



### Merging
Volume_total <- merge(Volume_taxation, Volume_unemployment, by = "Date") 


mi_overview<-Volume_total %>%
  dplyr::select(Date,taxationVolume,unemploymentVolume)%>%
  reshape2::melt( id.var='Date') %>%
  ggplot() + 
  geom_line( aes(x = Date, y = value, col=variable,alpha=0.3)) +
  scale_x_date( date_labels = "%Y",date_breaks = 'year') +
  xlab('Date') +
  ylab("coverage as % of total article count")+ 
  theme_bw()+
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"),
        legend.title = element_text(size=14), #change legend title font size
        legend.text = element_text(size=12))+
  scale_color_discrete(labels = c("taxation", "unemployment"), name = "Article phrase")+
  ggtitle("news coverage, Malawi")




news_benchmark<- ggarrange(ww_overview,mi_overview,
                      ncol = 2, nrow = 1,common.legend = TRUE, legend="bottom")

plot_media<- ggarrange(news_pre_rel1,news_pre_rel2,news_benchmark,
                          ncol = 1, nrow = 3)
plot_media

ggsave("2_output/fig_A18_media.pdf", width = 12, height = 12, units = "in")


```

## Figure A19: Malawi: Usage Frequency of Various Media
```{r , echo=TRUE,warning=FALSE}
# load afro-barometer 8 covariates
afro_8 <- read.spss("1_input/afrobarometer_release-dataset_merge-34ctry_r8_en_2023-03-01.sav", to.data.frame=TRUE)  %>%
   dplyr::filter(COUNTRY=="Malawi") %>%
  # select variables about news consumption
  dplyr::select(COUNTRY,RESPNO,Q55A,Q55B,Q55C,Q55D,Q55E)  %>%
  dplyr::rename(radio = Q55A,
                television = Q55B,
                newspapers = Q55C,
                internet = Q55D,
                social_media  = Q55E) 

# Now let us talk about the media and how you get information about politics and other issues. 55. How often do you get news from the following sources? 

long_afro_8 <- afro_8 %>%
  pivot_longer(cols = radio:social_media, names_to = "Media", values_to = "Frequency") %>%
  drop_na(Frequency) # Optionally drop NA values if any

# Calculate percentages
percent_data <- long_afro_8 %>%
  dplyr::group_by(Media, Frequency) %>%
  dplyr::summarise(Count = dplyr::n(), .groups = 'drop') %>%
  dplyr::group_by(Media) %>%
  dplyr::mutate(Percentage = Count / sum(Count) * 100)

fig_malawi_media<-ggplot(percent_data, aes(x = Frequency, y = Percentage, fill = Frequency)) +
  geom_bar(stat = "identity") +
  facet_wrap(~ Media, scales = "free_x") +
  theme_bw(base_size = 18) +
  labs(x = "", y = "Percentage", title = "Malawi: Usage Frequency of Various Media") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.position="none")
fig_malawi_media

ggsave("2_output/fig_A19_malawi_media.pdf", width = 14, height = 9, units = "in")

```


Figure A20: EU Speeches: Relative importance of topics between 2007 to 2015

```{r EU,warning=FALSE}
#######################################
# EU speeches
#######################################

load("1_input/euspeech.korpus.RData")


# Applying the Policy-Agendas und Laver-Garry dictionaries
load("1_input/policy_agendas_english.RData")
policyagendas.lexicon <- dictionary(dictLexic2Topics)


dfm.eu.pa <- dfm(korpus.euspeech, dfm_group = "country", dictionary = policyagendas.lexicon)



# Adding new categories with example keywords
dictLexic2Topics$disaster_preparedness <- c("preparedness","prevention","disaster preparedness","disaster prevention","emergency preparedness", "emergency plan", "disaster plan", "evacuation plan", "crisis management", "early warning system", "emergency kit", "emergency supplies", "disaster risk reduction", "resilience building", "community preparedness", "disaster preparedness training", "hazard mitigation", "emergency response plan", "preparedness exercises", "business continuity planning", "disaster resilience", "emergency shelter", "food storage", "water purification", "first aid training", "public health preparedness", "disaster simulation", "evacuation route", "emergency communication system")
dictLexic2Topics$disaster_relief <- c("relief","relief efforts", "aid distribution", "disaster recovery", "humanitarian aid", "emergency aid", "recovery operations", "disaster assistance", "rehabilitation efforts", "relief supplies", "medical relief", "search and rescue", "emergency medical services", "food aid", "shelter provision", "post-disaster recovery", "psychosocial support", "reconstruction assistance", "debris removal", "donation management", "volunteer coordination", "crisis counseling", "financial assistance", "temporary housing", "disaster victim identification", "infrastructure repair")

policyagendas.lexicon <- dictionary(dictLexic2Topics)

# Ensure your corpus has the necessary metadata, including the date
# This step is skipped here as it's assumed your corpus is already prepared

# Create a Document-Feature Matrix (DFM) using the dictionary for topic categorization
dfm.eu.pa <- dfm(korpus.euspeech) %>%
  dfm_lookup(dictionary = policyagendas.lexicon)


# Assuming 'Jahr' is the column in 'korpus.euspeech.stats' with the year information
# and it correctly aligns with the documents in your DFM
docvars(dfm.eu.pa, "year") <- korpus.euspeech.stats$Jahr


# Check for any NA values in year assignment
sum(is.na(docvars(dfm.eu.pa, "year")))

dfm_by_year <- dfm_group(dfm.eu.pa, groups = docvars(dfm.eu.pa, "year"))

# Convert counts to proportions within each year
dfm_prop <- dfm_weight(dfm_by_year, scheme = "prop")


# Convert DFM to a dense matrix, then to a data frame
dfm_matrix <- as.matrix(dfm_prop)
df <- data.frame(year = rownames(dfm_matrix), dfm_matrix)

# Use pivot_longer to convert wide data to long format
df_long <- pivot_longer(df, cols = -year, names_to = "topic", values_to = "frequency")

# Note: If the year column is not character, convert it to avoid issues in ggplot
df_long$year <- as.character(df_long$year)

target <- c("disaster_relief", "disaster_preparedness")

prep_relief<-df_long %>%
  filter(frequency<0.0025)%>%
  # filter( topic %in% target)  %>% # equivalently, dat %>% filter(name %in% target)
  ggplot( aes(x = year, y = frequency, color = topic, group = topic)) +
  geom_line() +
  theme_minimal() +
  labs(title = "Prepardness and Relief + similar importance",
       x = "Year",
       y = "Relative Importance",
       color = "Topic") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position="bottom")

all<-df_long %>%
  filter(frequency>0.05)%>%
  ggplot( aes(x = year, y = frequency, color = topic, group = topic)) +
  geom_line() +
  theme_minimal() +
  labs(title = "Top 5 topics",
       x = "Year",
       y = "Relative Importance",
       color = "Topic") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  theme(legend.position="bottom")


fig_eu_speech<- ggarrange(prep_relief,all,
                          ncol = 2, nrow = 1)
fig_eu_speech

ggsave("2_output/fig_A20_eu_speech.pdf", width = 12, height = 5, units = "in")

```
